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ABSTRACT 

Using the 330 hr ESO-MPG 870-^m survey of the Extended Chandra Deep Field South 
(ECDF-S) obtained with the Large Apex BOlometer CAmera (LABOCA) on the Atacama 
Pathfinder Experiment (APEX), we have carried out a stacking analysis at submillimeter 
(submm) wavelengths of a sample of 8266 near-infra-red (near-IR) selected (if vcga < 20) galaxies, 
including 893 BzK galaxies, 1253 extremely red objects (EROs) and 737 distant red galaxies 
(DRGs), selected from the Multi- wavelength Survey by Yale-Chile (MUSYC). We measure aver- 
age 870-/im fluxes of 0.20 ± 0.01 mJy (20.0a), 0.45 ± 0.04 mJy (11.3cr), 0.42 ± 0.03 mJy (U.Oa), 
and 0.41 ± 0.04 mJy (10.3a) for the if voga < 20, BzK, ERO and DRG samples, respectively. 
For the BzK, ERO and DRG subsamples, which overlap to some degree and are like to be at 
z ~ 1 — 2, this implies an average far-IR luminosity of ~ 2 — 6 x 10 11 L Q and star formation 
rate of ~ 40 — 100 M© . Splitting the BzK galaxies up into star-forming (sBzK) and passive 
(pBzK) galaxies, the former is significantly detected (0.48 ± 0.04 mJy, 12. Oct) while the latter is 
only marginally detected (0.27 ± 0.10 mJy, 2.7a), thus confirming that the sBzK/pBzK criteria 
do isolate obscured, star forming and truly passive galaxies. The K vcsa < 20 galaxies are found to 
contribute with 6.6 ± 0.3 Jy deg -2 (~ 15%) to the 870-/mi extragalactic background light (EBL) . 
sBzK and pBzK galaxies contribute 1.7 ± 0.2 Jy deg" 2 (-4%) and 0.2 ± 0.1 Jy deg" 2 (< 0.5%) 
to the EBL. We present the first delineation of the average submm signal from fC voga < 20 se- 
lected galaxies and their contribution to the submm EBL as a function of (photometric) redshift, 
and find a decline in the average submm signal (and therefore IR luminosity and star formation 
rate) by a factor ^2 — 3 from z ~ 2 to z ~ 0. This is in line with a cosmic star formation 
history in which the star formation activity in galaxies increases significantly at z > 1. A lin- 
ear correlation between the average 24-/zm and 870-/zm flux densities is found for K VBga < 20 
galaxies with 24-/itm fluxes < 350 fiJy (corresponding to L IR ~ 1.5 x 10 12 L az~2), while at 
higher 24-/mi fluxes there is no correlation. This behaviour suggests that star formation, and 
not Active Galactic Nuclei (AGN), is in general responsible for the bulk of the mid-IR emission 
of Li B £ 1-5 x 10 12 Lq systems, while in more luminous systems the AGN makes a significant 
contribution to the 24-^jm emission. By mapping the stacked 870-/im signal across the B — z vs. 
z — K diagram we have confirmed the ability of the sBzK-selection criterion to select starforming 
galaxies at z > 1, although our analysis suggest that the subset of sBzK galaxies which are also 
EROs are responsible for > 80 % of the submm emission from the entire sBzK population. 

Subject headings: cosmology: observations — galaxies: evolution — galaxies: formation — galaxies: 
high-redshift 



1. Introduction 

Extragalactic blank-field submm surveys have 
been carried out since the advent of SCUBA (Hol- 
land 1999) more than a decade ago, and have 
provided us with a unique view of intense, dust- 
cloaked star formation events at high redshifts 
(e.g. Blain et al. 2002). Yet such observations have 
so far only pin-pointed the most luminous high-z 
galaxies, due to the limitations in sensitivity and 
resolution imposed by preset-day (sub)mm facil- 
ities. It is now well-established that the bright 
(^ 5 mJy at 850-^m) submm sources uncovered by 
these surveys primarily reside in the redshift range 
z ~ 1.5 - 3.5 (Chapman et al. 2003, 2005), and 
account for ^ 20 — 30 percent of the extragalac- 
tic background light (EBL) at 850-/xm (Bargcr 
et al. 1998, 1999; Hughes et al. 1998; Coppin 
et al. 2006). Surveys that make use of galaxy 
clusters' gravitational amplification of the back- 
ground source plane have uncovered a number of 
faint (Ss50nm 2 mJy) sources and resolved up to 
80 per-cent of the background (Smail et al. 1997, 
2002a; Blain et al. 1999; Cowie et al. 2002; Chap- 
man et al. 2002; Knudscn et al. 2008). However, 
we know little about the nature and redshift dis- 
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tribution (Smail et al. 1997, 2002a) of the popula- 
tion below ~ 5 mJy due in part to the difficulty of 
identifying counterparts in the radio. 

The recent advent of large format near-IR cam- 
eras have revealed populations of moderately star 
forming galaxies at z ~ 1 — 3 that are more numer- 
ous than the (sub)mm selected systems, and more 
representative of the bulk population at these 
epochs (e.g. Cimatti et al. 2002; Lawrence et al. 
2007). The rest- frame near-IR is arguably the best 
wavelength range to undertake such surveys at - 
as, in comparison with UV and optical surveys, it 
is less sensitive to the effects of age and dust on the 
stellar population, and thus more closely provides 
a selection based on stellar mass. 

Well known examples of near-IR colour-selected 
galaxies are the populations of extremely red ob- 
jects (EROs - Elston, Rieke & Rieke 1988, Mc- 
Carthy, Persson & West 1992; Hu & Ridgway 
1994), distant red galaxies (DRGs - Franx et 
al. 2003; van Dokkum et al. 2003), and the so- 
called BzK galaxies (Daddi et al. 2004). While 
these populations are selected according to differ- 
ent optical/near-IR colour criteria applied, which 
pick out systems at different, but overlapping, 
redshift ranges. These colour criteria are often 
designed to straddle the 4000 A break (includ- 
ing the Balmer break at 3646 A), characteristic 
of evolved, metal-enriched galaxies that are old 
enough that OB-stars do not dominate the light. 
The same colour criteria, however, will also select 
dusty, starforming galaxies at virtually any red- 
shift. Thus, near-IR colour-selected galaxy popu- 
lations are typically a mix of actively star forming 
galaxies and old, evolved systems, which means 
additional colour criteria and/or spectral analysis 
has to be applied in order to separate the two. 
Clearly, submm observations offer a unique way 
of distinguishing between starforming and passive 
near-IR galaxies. 

At present, however, the bulk of near-IR colour- 
selected galaxies arc too faint for individual de- 
tection by large format (sub)mm surveys, and for 
the moment, therefore, the only way forward is 
to study their average (sub)mm/far-IR proper- 
ties by means of stacking techniques. A hand- 
ful of such studies have been carried out to date, 
characterizing the average submm signal of near- 
IR selected galaxies and their contribution to the 
extragalactic background light (EBL) at submm 
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wavelengths (Webb et al. 2004; Daddi et al. 2005; 
Knudsen et al. 2005; Takagi et al. 2007). Yet most 
of these stacking analyses have been of relatively 
small samples of galaxies, and as a consequence 
have had to averaging their submm properties over 
the entire redshift range from which they are se- 
lected (which is often substantial, z ~ 1 — 3). A 
robust delineation of the submm signal of near- 
IR selected galaxies as a function of redshifts has 
therefore been lacking, and as a consequence we do 
not know how the dust-enshrouded starformation 
in these galaxies evolve with cosmic epoch. 

The Extended Chandra Deep Field South 
(ECDF-S), a 30' x 30' region centered on the 
smaller GOODS-S/CDF-S field (Giavalisco et 
al. 2004), is one of the most intensively stud- 
ied extragalactic fields in the southern sky. In 
addition to X-ray observations with Chandra 
(Alexander et al. 2003; Lehmer et al. 2005; 
Luo et al. 2008), the EDCF-S has been tar- 
geted in a large number of optical and near- 
IR filter passbands from the ground as part of 
COMBO-17 (Wolf et al. 2001) and MUSYC (Ga- 
wiser et al. 2003), and with HST/ACS as part 
of GEMS (Rix et al. 2004). Furthermore, deep 
mid-IR imaging has been provided by the Spitzer 
IRAC/MUSYC Public Legacy in ECDF-S (SIM- 
PLE, Damen et al. 2009), and the Spitzer /MIPS 
Far-Infrared Deep Extragalactic Legacy Survey 
(FIDEL - Dickinson et al. in prep. See also 
http : //ssc . spitzer . caltech.edu/legacy/abs 
/dickinson2 . html). 

To study the submm properties of the sources 
in the ECDF-S, we have undertaken a large ESO- 
MPG survey (Coppin et al. 2009; Weifi et al. 2009) 
using the LABOCA 870-/im camera (Siringo et al. 
2009) mounted on the Atacama Pathfinder Exper- 
iment (APEX) and combined the data with the al- 
ready existing multi-wavelength data available for 
this field. 

Throughout this paper we adopt a flat cosmol- 
ogy with M = 0.27, n A = 0.73, and h = 0.71 
(Spergel et al. 2003). 

2. The submm data 

Observations were carried out using the 295 
horn-bolometer array LABOCA on APEX (Siringo 
et al. (2009), and are discussed in detail in Weifi et 
al. (2009)). The bolometers are AC-biased, op- 



erated in total power mode, and distributed in 
a hexagonal configuration over the 11.4' field 
of view. The center frequency of LABOCA 
is 345 GHz and its passband has a FWHM of 
~ 60 GHz. The measured angular resolution is 
19.2" (FWHM). The observations were carried out 
between May 2007 and November 2008 in mostly 
excellent weather conditions (PWV typically 0.5 
mm corresponding to a zenith opacity of 0.2 at the 
observing wavelength) . The mapping pattern was 
chosen to give a uniform coverage across a 30' x 30' 
area centered at Ra: 03' l 32 m 29 ;i Dec -27°48'47" 
(J2000). 

Mapping was performed alternating rectangu- 
lar on-the-fly scans with raster of spiral patterns. 
For the latter mode the telescope traces in two 
scans spirals with radii between 2' and 4' at 16 and 
9 positions (the raster) spaced by 600" in azimuth 
and elevation. The scanning speed was typically 
between 2-3 arcmin per second. 

Calibration was achieved through observations 
of Mars, Uranus and Neptune as well as secondary 
calibrators and was found to be accurate within 
8.5%. The atmospheric attenuation was deter- 
mined via skydips every ~ 2 hours as well as 
from independent data from the APEX radiome- 
ter which measures the line of sight water va- 
por column every minute. Focus settings were 
typically determined once per night and checked 
during sunrise. Pointing was checked on nearby 
quasars PMNJ0457-2324, PMNJ0106-4034 and 
PMNJ0403-3605 and found to be stable within 
3". 

The data was reduced using the BoA reduc- 
tion package (Schuller et al. in prep.). Individual 
maps were co-added (noise- weighted) and the final 
map was beam-smoothed, resulting in a spatial 
resolution of 27" (FWHM). The total on-source 
observing time in the data used for this analy- 
sis is 200 hours (330 hr including overheads) and 
the average rms across the entire 30' x 30' field 
is 1.2 mjy beam -1 , making it the largest contigu- 
ous (sub)mm survey ever undertaken to this depth 
(c.f. Coppin et al. 2006; Bertoldi et al. 2007). 

3. Near-IR selected galaxies 

We use the Wide MUSYC public data release of 
UBVRIzJHK catalogues in the ECDF-S to con- 
struct samples of near-IR selected galaxies (Taylor 
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et al. 20080 The MUSYC survey covers the cen- 
tral 30' x 30' of the ECDF-S and reaches 5-cr point 
source sensitivities in the near-IR of J AB = 22.7 
and 

4. Near-IR selected galaxies 

We use the Wide MUSYC public data release of 
UBVRIzJHK catalogues in the ECDF-S to con- 
struct samples of near-IR selected galaxies (Taylor 
et al. 2008 jE The MUSYC survey covers the cen- 
tral 30' x 30' of the ECDF-S and reaches 5-cr point 
source sensitivities in the near-IR of J AB = 22.7 
and K AB = 22.0, respectively (see Gawiser et al. 
(2006) for a detailed description of the MUSYC 
survey) . 

We defined our sample as sources with if voga < 
20 (corresponding to K AB < 21.9, i.e. we have 
assumed a VEGA-AB offset of 1.9 in the K- 
band). This magnitude cut-off was chosen since 
the MUSYC survey is close to 100% complete at 
this depth, and since other studies have adopted 
the same magnitude limit, thus facilitating a di- 
rect comparison. Contamination by stars was ac- 
counted for by removing objects lying on the stel- 
lar loci in the (z — K) AB vs. (B — z) AB or (J — K) AB 
vs. (R — K) AB diagrams (hereafter referred to as 
BzK and RJK diagrams - Fig. [I]). As shown in 
Fig. QJ, and b, a non-negligible fraction of sources 
that appear as normal galaxies in the BzK diagram 
end up in the stellar region in the RJK diagram 
and vice versa. These sources were all poorly fit 
(xV^dot > 100) by galaxy SED templates, thus 
strongly suggesting that they are in fact stars and 
illustrating the need of removing both stellar loci 
from the sample. This left us with a total sample 
of 8266 iCcga < 20 selected galaxies. From this 
sample, subsets of BzK galaxies, extremely red 
objects (EROs) and distant red galaxies (DRGs) 
were extracted as described below. 

Applying the star-forming and passive BzK cri- 
teria, namely BzK = (z — K) AB — (B — z) AB > 
—0.2 for sBzK galaxies, and BzK < —0.2 and 
(z — K) AB > 2.5 for pBzK galaxies (Daddi et al. 
2004), to our A' voga < 20 catalogue defined above, 
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we obtained samples of 744 sBzK and 149 pBzK 
galaxiedE- Sources which were undetected at the 
1 a level in both B and z were not included. Fur- 
thermore, sources formally belonging to the sBzK- 
region but which were undetected in B, thus hav- 
ing B — z lower limits only, were discarded, as 
were sources in the pBzK region with lower lim- 
its in z — K. A total of 1253 EROs were selected 
from our -ftT voga < 20 sample using the standard 
criterion, i.e. (R — K) AB > 3.35 (equivalent to 
(R - if) voga £ 5 - Elston, Rieke & Rieke 1988; 
McCarthy, Persson & West 1992; Hu & Ridgeway 
1994). Finally, the selection of DRGs was done by 
applying the colour cut (J — K) AB > 1.32 (Franx 
et al. 2003). In total 737 DRGs with isT vega < 20 
were selected in this way. 

4.1. Overlap between BzK, ERO, and 
DRG samples 

In this paper we aim to determine the contri- 
bution to the submm background from if voga < 20 
galaxies, as well as the sub-samples of BzKs, EROs 
and DRGs and their joint contributions, and it 
is therefore important to determine the degree of 
overlap between these populations. The overlaps 
in terms of percentages are given in Table [T] 

From Fig. |2] it is seen that BzK, ERO and 
DRG galaxies only start to contribute significantly 
(> 1%) to the full sample for K AB £ 20. Even 
so, more than half of the full sample does not 
fall within the BzK/ERO/DRG classifications at 
these faint flux levels. Of the full if voga < 20 sam- 
ple, 6269 sources (corresponding to ~ 76%) do not 
classify as BzKs, EROs or DRGs. 

The locations of the various samples in the BzK 
and RJK diagrams are shown in Fig. [1] The DRGs 
are seen to be spread out across the BzK diagram, 
while the EROs lie in a much more well defined re- 
gion of the BzK diagram. EROs and DRGs make 
up about 30% and 36% of the sBzKs, respectively. 
Similarly, we find that EROs and DRGs constitute 
98% and 44% of the pBzK sample, respectively. 
Clearly, pBzK galaxies are much more often se- 
lected as EROs than is the case for sBzKs, while 
the occurrence of a pBzK or sBzK galaxy being 
classified as a DRG is about the same. Turning 
to the RJK-diagram we see that the sBzK galax- 
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Table 1 

The total number of sources in each sample (given in parentheses in the first column of each line), 
followed by the percentages contributed by the other samples. the latter are also given as 

absolute numbers in parentheses. 





sBzK 


pBzK 


ERO 


DRG 


sBzK+pBzK 


ERO+DRG 


if voga < 20 (8266) 


9%(744) 


1.8%(149) 


15.296(1253) 


8.9%(737) 


10.8%(893) 


19.696(1620) 


sBzK (744) 


100%(744) 


0%(0) 


30.496(226) 


36.2%(269) 


83.396(744) 


49.796(370) 


pBzK (149) 


0%(0) 


100%(149) 


98.0%(146) 


43.6%(65) 


16.796(149) 


9896(146) 


ERO (1253) 


18.0%(226) 


11. 7% (146) 


10096(1253) 


29.596(370) 


29.796(372) 


10096(1253) 


DRG (737) 


36.596(269) 


8.8%(65) 


50.2%(370) 


100%(737) 


45.396(334) 


100% (737) 


sBzK+pBzK (893) 


83.396(744) 


16.796(149) 


41.796(372) 


37.496(334) 


100% (893) 


57.896(516) 


ERO+DRG (1620) 


22.896(370) 


9.0%(146) 


63.0%(1253) 


37.0%(737) 


31.9%(516) 


10096(1620) 



ies are much more spread out than the pBzKs. 
About 18% and 37% of EROs and DRGs, respec- 
tively, are made up of sBzK, while pBzKs make 
up about 12% and 9% of the two populations. 

The overlaps between BzK, ERO and DRG 
galaxies have been discussed in detail in other 
studies (e.g. Reddy et al. 2005; Grazian et al. 2007; 
Takagi et al. 2007; Lane et al. 2008). The most 
statistically significant study was carried out by 
Lane et al. (2009) who used the UKIRT Infrared 
Deep Sky Survey (UKIDSS) Ultra Deep Survey 
Early Data Release (UDS EDR) to study large 
samples of BzK, ERO and DRG galaxies. For sam- 
ples selected down to K AB = 21.2, which is close 
to our magnitude cut-off, they found sBzK:ERO 
and pBzK:ERO ratios of 32% and 95%, respec- 
tively, i.e. in excellent agreement with our values. 
They also find that about 30% of DRGs are sBzK, 
again in good agreement with our findings (see 
also Reddy et al. 2005). 

4.2. Redshift distributions 

The sample was correlated against publically 
available spectroscopic redshift surveys (Szokoly 
et al. 2004; Vanzella et al. 2005, 2006, 2008; 
Popesso et al. 2008; Kriek et al. 2008). Using 
only the most reliable spectroscopic redshifts from 
these surveys we extracted 2341 redshifts, of which 
the majority lie within the CDF-S region. A total 
of 546 galaxies from our sample were matched to a 
spectroscopic redshift. Of these 28 were sBzK, 21 
were ERO, and 10 were DRG galaxies. The fact 
that no pBzK galaxies were identified with a spec- 
troscopic redshift is not too surprising since these 
are in all likelihood old, evolved galaxies with opti- 
cal spectra devoid of emission features, thus mak- 



ing it difficult to obtain robust spectroscopic red- 
shifts (cf. Kriek et al. 2006; Cimatti et al. 2008). 
Sources with spectroscopic redshifts were used as 
a test sample to optimize input parameters for the 
photometric redshift code EAZY (Brammer, van 
Dokkum & Coppi 2008). The code works by fitting 
non-negative linear combinations of galaxy spec- 
tra to the observed spectral energy distributions 
(SEDs), which in our case consisted of the nine 
MUSYC filter bands. 

The resulting photometric redshifts are com- 
pared against their spectroscopic counterparts in 
Fig. [3] The normalized median absolute deviation 
of Az = z phot — z 3pcc (see Brammer, van Dokkum 
& Coppi 2008) is ~ 0.037 for z < 1.5 and ~ 0.079 
for z > 1.5. Significant outliers, which we define 
to be sources with \Az\/ (1 + z apo[ .) 5 times greater 
than the median, make up ~ 9% of the total sam- 
ple. These numbers are consistent with the typical 
performances of photometric redshift codes (e.g. 
Bolzonella, Miralles & Pello 2000; Quadri et al. 
2007; Brammer van Dokkum & Coppi 2008). 

Adopting the parameters for the test sample, 
photometric redshifts were derived for the remain- 
der of the sample without spectroscopic redshifts. 
The redshift distributions of the full K voga < 20 
sample as well as the BzK, ERO, and DRG sam- 
ples are shown in Fig. Q] where we also com- 
pare with two other photometric redshift distri- 
butions obtained from: 1) the MUSYC UBVRI 
and deep JHK imaging of the three 10' x 10' fields 
HDFS1/S2, MUSYC 1030 and 1255 (Quadri et 
al. 2007 - hereafter Q07), and 2) deep BVRi'z'JK 
imaging of 1113arcmin 2 in the Ultra-Deep Survey 
(UDS) portion of the United Kingdom Infrared 
Telescope Deep Sky Survey (UKIDSS) (Dunne et 
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Fig. 1. — Colour-magnitude diagrams illustrating the selection criteria for the BzK, ERO and DRG samples. 
Left: The BzK-diagram (D04) showing the sBzK and pBzK regions and the location of sources selected 
according to the ERO- and DRG-criteria. Right: RJK diagram showing the ERO and DRG regions and the 
location of the sources fulfilling the sBzK and pBzK colour criteria. Sources below the dotted lines in either 
diagram were deemed to be stars and therefore discarded from the analysis. Note this includes a number of 
sources formally selected as either BzK, ERO or DRG galaxies (see §|4j. The diagrams show the significant 
overlap between the different near-IR selected populations. 



al. 2008 - hereafter D08). 

The BzK selection criteria are defined to se- 
lect galaxies in the redshifts range 1.4 < z < 2.5 
(Daddi et al. 2004), yet both the sBzK and pBzK 
redshift distributions extend below and above this 
range. Of the sBzK galaxies, 66 % lie in the range 
1.4 < z < 2.5, while 22% are at z < 1.4 and 12% 
at z > 2.5. For the pBzK galaxies, the correspond- 
ing percentages are 60%, 31%, and 9%, respec- 
tively. Both Q07 and D08 find similar fractions 
for their samples of if vega < 20 galaxies, suggest- 
ing that the BzK criteria select galaxies across a 
somewhat broader redshift range (1 £ z £ 3.5). 
The redshift distribution of EROs is seen to peak 
strongly at z ~ 1.1 with a tail extending to 
z ~ 3.5. This is in line with photometric and spec- 
troscopic surveys which have shown that the red- 
shift distribution of EROs peaks around z ~ 1.2 
(Cimatti et al. 2003; Yan et al. 2004). The red- 
shift distribution of if voga < 20 EROs derived by 
D08 peaks at slightly higher redshifts (z ~ 1.4), 



but overall appears similar to ours. The DRG 
distribution shows prominent peaks at z ~ 1.2 
and z ~ 2, with the former being the most dom- 
inant. A similar bimodality is also apparent in 
the X vcga < 20 DRG sample by Q07, although the 
dominant peak in their distribution lies at z ~ 2. 
The distribution by D08 broadly resembles ours, 
with a prominent peak at z ~ 1.2 followed by a 
significant high-z tail. Overall, therefore, our red- 
shift distributions are consistent with those of Q07 
and D08, given the uncertainties associated with 
photometric redshift derivation, and the effects of 
field-to-field variations. 

5. 870-^m Stacking 

In order to estimate the average 870-/im fluxes 
for the above excised catalogues of near-IR galax- 
ies we stack the 870-/im flux values at their near- 
IR positions in the LABOCA map. 
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Fig. 2. — Top: Percentages, i.e. fraction of the 
full if voga < 20 sample, of the BzK, ERO and 
DRG samples, as well as the non-BzK/ERO/DRG 
sample (dashed line) as a function of their K A b 
magnitude. Bottom: The number distribution 
of sources as a function of their K AB magnitude. 
The distributions for the full sample as well as the 
BzK, ERO and DRG samples are shown in order 
to illustrate the overlap between the different pop- 
ulations. While the BzK, ERO and DRG samples 
start making up a non-negligible fraction of the 
parent sample for K AB > 20, even at the faintest 
magnitudes more than half of the sources from the 
parent sample do not fall within any of these three 
classifications. 

5.1. Submm-bright near-IR selected galax- 
ies 

First, however, we need to identify any galax- 
ies which are associated with robust LABOCA 
sources, which we take to mean sources detected 
at > 3.7-ct significance (see Weifi et a!. 2009 for 
details). To this end we adopt a search radius 
of 12.8"around each LABOCA source, which cor- 
responds to the 95 per-cent confidence search ra- 
dius given the FWHM=19" beam. For a galaxy 
to qualify as a near-IR counterpart to a LABOCA 
source, we furthermore require that z > 0.8. If 
more than one galaxy meets these criteria for a 
given LABOCA source, we adopt the one closest 
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Fig. 3. — Photometric versus spectroscopic red- 
shifts for the sources in the full i^ veg a < 20 sam- 
ple, and the sBzK/pBzK, ERO and DRG sub- 
samples which have spectroscopic redshifts (see 
§ I4.2j) . The normalized median absolute deviation 
of Az — z phot — z spoc (see Brammer, van Dokkum 
& Coppi 2008) is ~ 0.037 for z < 1.5 and ~ 0.079 
for z > 1.5. 

to the submm source. In this manner we find 54 
submm-near-IR associations from the K vcsa < 20 
sample. Of these, 17/2 sources are classified as 
sBzK/pBzK galaxies and 28 as EROs (of which 
16 are also DRGs). 11 of the 17 sBzK galaxies 
are also EROs and of those, 9 are DRGs. Both 
of pBzK galaxies are EROs, and also DRGs. As 
a safeguard against contamination, these submm- 
bright near-IR sources were removed from the 
stacking analysis, although their contribution was 
included (in a variance-weighted fashion) in the 
final tally of average submm flux (Table [5]). 

5.2. Stacking and deblending technique 

Next, we proceeded to perform a stacking anal- 
ysis of the remaining submm-undetected near-IR 
selected galaxies. Due to the slightly varying noise 
across the map, the average 870- /im flux ((SWo^m)) 
and noise ((<7 87 o Mm )) values were calculated as the 
variance-weighted mean, i.e. 



and 



(Cs70^m) = 



i 



(1) 

(2) 
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Fig. 4. — Top to bottom: redshift distributions of 
the X vcga < 20, sBzK, pBzK, ERO and DRG sam- 
ples, respectively (dark grey histograms). Spectro- 
scopic redshifts are shown as light grey histograms. 
Also shown are the photometric redshift distribu- 
tions of the JT voga < 20 sBzK, pBzK and DRG sam- 
ples by Q07, as well as of the -ftT voga < 20 selected 
samples by D08. The Q07 and D08 distributions 
have been scaled to match the peaks of our distri- 
butions. Given the uncertainties, our if voga < 20 
BzK, ERO and DRG redshift distributions are in 
good agreement with those of Q07 and D08. 

where Si and cr, are the 870-/im flux and r.m.s. 
noise pixel values at the near-IR position of the 
i th source in the stack, respectively. To avoid 
the stack being contaminated by robust 870-/zm 
sources, the stacking was performed on a 'residual' 
version of the LABOCA map in which all of the 
55 870-^m sources uncovered (Weifi et al. 2009) 
had been subtracted using a scaled beam profile. 

An important aspect of any stacking analysis 
performed on maps with coarse angular resolution, 
is the issue of deblending of sources that lie within 
a single resolution element. For example, if a BzK 
galaxy has a neighbour, A, within a LABOCA 



beam, we have to calculate the 870-^m flux contri- 
bution from A at the position of the BzK galaxy. 
Of the full -ftT voga < 20 sample, we find this to be 
the case for 5985 sources (i.e. 72 per-cent of the 
sample). We correct for the blending of sources 
by assuming Gaussian sources with FWHMs equal 
to the LABOCA beam. This method is similar to 
the one adopted by Webb et al. (2004), although 
they (and subsequent submm stacking studies) did 
not take into account the effects from neighbour's 
neighbours. An illustrative example of the lat- 
ter is given Fig. [5] where A itself has a neighbour, 
B, within a LABOCA beam, that is not within 
the LABOCA beam as measured from the BzK 
galaxy's position. In this case we have to calcu- 
late B's contribution to A, in order to correctly 
calculate A's contribution to the BzK galaxy, and 
the system of equations to be solved is therefore 

/b,k = I Bz K+/ A e- r -K,A/(2- 2 ) ( 3 ) 

h = I A + I B ^e-<^/^+I B e- r l^^) 
f B = I B +I A e- r l-/^ (5) 

where I and / are the measured and deblended 
fluxes at the relevant positions, respectively, and 
r are the distances between the sources. In order 
to estimate the error one makes by only deblend- 
ing the fluxes from neighbours within a beam, 
not taking into account neighbours' neighbours, 
we ran the stacking analysis under both scenarios. 
We find that the deblending scheme by Webb et 
al. (2004) can overestimate the fluxes by ~ 10% 
compared to the scheme described in this pa- 
per, but more typically the error is at the ~ 5 % 
level. We emphasize that one has to not only 
dcblcnd neighbours within a certain galaxy pop- 
ulation, but also across populations. Therefore, 
the deblending analysis was carried out on the full 
-?Ce g a < 20 sample. In doing so we have ignored 
sources with K VBga > 20, and they have therefore 
not been included in the deblending analysis. In 
the following section, however, we show that these 
fainter sources do indeed make a contribution to 
the submm signal, which must be subtracted from 
the stacked fluxes. 

5.3. Stacking the full samples 

For each source in the -ftT VBga < 20 sample the 
(deblended) signal and noise values at its pixel po- 
sition in the LABOCA map were recorded, and 



S 



BzK 



B 



ysis. As expected the distributions of stacked sig- 
nals obtained from these simulations were Gaus- 
sians centered on zero. We found that the mea- 
sured 870-/mi signals occurred in < 0.05% of the 
simulation runs. Roughly the same percentage is 
found if we restrict our Monte Carlo analysis to the 
sBzK, ERO and DRG samples, and clearly testify 
to the significance of the measured signals. 



Fig. 5. — An example of where deblending is nec- 
essary. In order to properly calculate the 870-/im 
flux coming from the BzK galaxy at its position we 
have to calculate A's contribution, which in turn 
is affected by B. In total, therefore, a linear set 
of three equations has to be solved (see eqs. 3-5). 
The large circles indicate the LABOCA beam. 



from those the stacked 870-/im flux density of the 
full sample was determined according to eqs. [T] 
- [2J The 870- fim signal and noise values corre- 
sponding to the BzK, ERO and DRG samples were 
extracted and their stacked 870-^m flux densities 
were derived in a similar manner. In a similar way, 
postage stamp images around each source were 
extracted and combined in a weighted fashion re- 
sulting in stacked submm-images of the if-selected 
samples (Fig. [B]). From the azimuthally-averaged 
radial profiles of the submm signal, it is clear that 
the baseline level is not zero, but in fact there is a 
residual signal amounting to 0.065 mJy and stem- 
ming from the population of K vcga > 20 galax- 
ies lying below the submm detection limit of the 
LABOCA map. In order to account for this effect, 
the final stacked 870-^m flux densities (Table [2]) 
had a constant signal of 0.065 mJy subtracted from 
them. As a comparison we also derived the median 
flux densities from the stacks and found agreement 
(to within 15 percent) with the weighted averages. 

To gauge the significance of our results, we ran 
a series of Monte Carlo simulations in which stack- 
ing analyses were carried out on 1000 versions of 
the -fCoga < 20 catalogue, each with randomized 
positions with respect to the original catalogue. 
Each source was assigned a random position by 
(randomly) choosing a radius (60" < r < 200") 
and an angle from its original position. By confin- 
ing the new positions within a certain distance of 
the original, we ensured that the noise properties 
were similar to those in the original stacking anal- 
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Fig. 7. — The distributions of the S/N values 
at the positions of the # Tega < 20 BzK, ERO, 
and DRG samples, shown as hashed histograms 
from top to bottom, respectively. The open his- 
tograms show the S/N distribution of the residual 
LABOCA map. The P KS values arc the likelihoods 
that the sample distributions are identical to the 
overall map distribution. 

In addition to the above stacking analysis, we 
compared the distributions of S/N values at the 
near-IR source positions with the overall S/N 
distribution of the residual LABOCA map (Fig. 
[7]). In all cases, a formal Kolmogorov-Smirnov 
(KS) statistic shows that the sample distributions 
are statistically different from the distribution of 
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Table 2 

The stacked 870-/im flux densities for the K^ ECA < 20 sample as well as the BzK, ERO, and DRG 
samples. The numbers of sources going into each stack are given in parentheses along with the 
significance of the stacked signals. surface densities and contribution to the 870-/1m ebl are also 

LISTED. 



Galaxy type 


(S870(Jm) a 


(S870(Jm) b 


<5N b 


(<5Is70nm> b 




[m 


Jy] 


[m 


Jy] 


[sq. arcmin -1 ] 


[Jy sq. dog -1 ] 


ifvcga < 20 


0.19 ± 0.01 


(8212, 19. Off) 


0.20 ± 0.01 


(8266, 20. Oct) 


9.18 ± 0.10 


6.61 ± 0.34 


(15.0 ± 5.2%) 


sBzK 


0.38 ± 0.04 


(727, 9.5<t) 


0.48 ± 0.04 


(744, 12. Oct) 


0.83 ± 0.03 


1.70 ± 0.24 


(3.9 ± 1.4%) 


pBzK 


0.21 ± 0.10 


(147, 2.1ct) 


0.27 ± 0.10 


(149, 2.7ct) 


0.17 ±0.01 


0.16 ± 0.11 


(0.4 ±0.3%) 


sBzK+pBzK 


0.35 ± 0.04 


(874, 8.8a) 


0.45 ± 0.04 


(893, 11.3<7) 


0.99 ± 0.03 


1.89 ± 0.24 


(4.3 ± 1.6%) 


ERO 


0.30 ±0.03 


(1225, 10. Oct) 


0.42 ± 0.03 


(1253, 14. Oct) 


1.39 ± 0.04 


1.85 ± 0.20 


(4.2 ± 1.5%) 


DRG 


0.32 ± 0.04 


(717, 8. Oct) 


0.41 ± 0.04 


(737, 10.3ct) 


0.82 ± 0.03 


1.09 ± 0.18 


(2.5 ±0.9%) 


BzK/ERO/DRG c 


0.31 ± 0.03 


(1961, 10.3ct) 


0.43 ± 0.03 


(1997, 14.3ct) 


2.22 ± 0.05 


3.67 ±0.29 


(8.4 ± 2.9%) 


non-BzK/ERO /DRG d 


0.16 ± 0.01 


(6251, 16. Oct) 


0.17 ± 0.01 


(6269, 17. Oct) 


6.97 ± 0.09 


4.26 ± 0.26 


(9.7 ± 3.3%) 








2 < 


1.4 








ifvcga < 20 


0.17 ± 0.01 


(7026, 17. Oct) 


0.17 ± 0.01 


(7055, 17. Oct) 


7.83 ± 0.09 


5.36 ± 0.29 


(11.7 ±4.1%) 


sBzK 


0.33 ± 0.09 


163, 3.7ct) 


0.42 ± 0.09 


(164, 4.7a) 


0.18 ± 0.01 


0.23 ± 0.14 


(0.5 ± 0.4%) 


pBzK 


-0.03 ± 0.17 (47, -0.18ct) 


-0.14 ± 0.17 (47, -0.8a) 


0.05 ± 0.01 




ERO 


0.24 ± 0.04 


(740, 6.0ct) 


0.26 ± 0.04 


(751, 6.5ct) 


0.83 ± 0.03 


0.87 ± 0.16 


(1.9 ± 0.7%) 


DRG 


0.26 ± 0.06 


(351, 4.3a) 


0.29 ± 0.06 


(373, 4.8a) 


0.40 ± 0.02 


0.47 ± 0.15 


(1.0 ± 0.4%) 








2 > 


1.4 








ifvcga < 20 


0.34 ± 0.03 


(1186, 11.3<r) 


0.39 ± 0.03 


(1211, 13. Oct) 


1.35 ± 0.04 


1.89 ± 0.20 


(4.1 ± 1.5%) 


sBzK 


0.40 ± 0.05 


(564, 8. Oct) 


0.55 ± 0.05 


(580, 11. Oct) 


0.64 ± 0.03 


1.28 ± 0.23 


(2.8 ± 1.1%) 


pBzK 


0.32 ± 0.12 


(100, 2.7a) 


0.41 ± 0.12 


(102, 3.4ct) 


0.11 ± 0.01 


0.17±0.16 


(0.4 ± 0.4%) 


ERO 


0.40 ± 0.05 


(485, 8. Oct) 


0.54 ± 0.05 


(502, 10.8ct) 


0.56 ± 0.02 


1.08 ± 0.22 


(2.4 ±0.1%) 


DRG 


0.39 ± 0.06 


(366, 6.5a) 


0.59 ± 0.06 


(380, 9.8ct) 


0.42 ± 0.02 


0.90 ± 0.23 


(1.9 ± 0.8%) 




24-/^m detected (S24^m > 27 /^Jy) 


sBzK 


0.55 ± 0.05 


(451, 11. Oct) 


0.70 ± 0.05 


(466, 14. Oct) 


0.52 ± 0.02 


1.30 ± 0.27 


(2.8 ± 1.1%) 


pBzK 


0.44 ± 0.16 


(52, 2.8a) 


0.54 ± 0.16 


(53, 3.4ct) 


0.06 ± 0.01 


0.11 ± 0.20 


(0.2 ± 0.4%) 


ERO 


0.54 ± 0.05 


(490, 10.8ct) 


0.68 ± 0.05 


(511, 11.6ct) 


0.57 ±0.03 


1.39 ± 0.28 


(3.0 ± 1.0%) 


DRG 


0.52 ± 0.06 


(379, 8.7a) 


0.73 ± 0.06 


(400, 12.2ct) 


0.44 ± 0.02 


1.16 ± 0.28 


(2.5 ± 1.0%) 








24- /^m faint 


(S 2 4 M m < 27AiJy) 








sBzK 


0.11 ± 0.07 


(276, 1.6a) 


0.14 ± 0.07 


(278, 2. Oct) 


0.31 ± 0.02 


0.16 ± 0.09 


(0.3 ± 0.2%) 


pBzK 


0.08 ± 0.12 


(95, 0.7ct) 


0.12 ± 0.12 


(96, I.Oct) 


0.11 ± 0.01 


0.05 ± 0.06 


(0.1 ±0.1%) 


ERO 


0.14 ±0.04 


(735, 3.5a) 


0.18 ± 0.04 


(742, 4.5a) 


0.82 ± 0.03 


0.53 ± 0.14 


(1.2 ±0.5%) 


DRG 


0.10 ±0.06 


(338, 1.7a) 


0.18 ± 0.06 


(344, 3.0ct) 


0.38 ± 0.02 


0.25 ± 0.11 


(0.5 ± 0.2%) 



Note. — a Submm-faint sources only. b Submm-bright sources included. c The combined BzK, ERO and DRG samples, where source 
overlap between the populations have been accounted for. 
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Fig. 6. — Top row: Postage stamps of the stacked signal-to-noise images, each 310" x 310" in size, of the 
if lcga < 20, BzK, ERO and DRG samples. The contours start at S/N = 3 and increase in steps of 2. The 
number of sources in each stack is given in parentheses. The angular resolution (FWHM=27") is shown as 
an insert in the left hand panel. Middle row: Radial profiles (azimuthally averaged) of the corresponding 
stacked S/N images (filled symbols). Gaussian fits to the profiles are indicated by the solid curve. Bottom 
row: Radial profiles (azimuthally averaged) of the corresponding stacked signal images (filled symbols). 
Gaussian fits to the profiles are indicated by the solid curve. Notice the non-zero baseline level (0.065 mJy) 
caused by the the K vcsa > 20 sources with submm fluxes below the detection threshold. This constant 
baseline level has been removed from the S/N images and profiles in the top and middle panels. 



the entire map which is symmetrical around zero 
(open histogram in Fig. [7|). Unlike the distribu- 
tions for the £T vcga < 20, sBzK, ERO and DRG 
samples, which are clearly biased towards positive 
S/N values, the pBzK distribution shows outliers 
at both positive and negative S/N values. We can 
therefore not confidently claim a detection even 
though this is formally suggest by the low KS 
probability. This is in line with the above anal- 
ysis which found that the stacked 870-/im signal 
from pBzKs was not significant. 



6. Discussion 

6.1. Stacked submm fluxes and star for- 
mation rates 

In §[5] we found that all of our if-selected sam- 
ples, except for the pBzKs, have significantly de- 
tected stacked 870-/im fluxes. The submm-faint 
pBzK galaxies are tentatively detected at the ~ 
2cr-level - as also indicated by the slight increase in 
the azimuthally averaged S/N radial profile in Fig. 
[7] How do these compare with previous submm 
stacking analyses of if-selected samples? 

We find average 870-/im flux densities of 0.38 ± 
0.04 mJy and 0.21 ± O.lOmJy for our samples of 
submm-faint sBzK and pBzK galaxies. In com- 
parison, Daddi et al. (2005) estimated an aver- 
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age 850-^m signal of ~ 0.82 mJy (corresponding 
to a 870-/xm signal of ~ 0.75 mJjo) for a sample 
of ~ 100 if voga < 20 submm- faint sBzK galax- 
ies within the SCUBA map of GOODS-N. Tak- 
agi et al. (2007) reported an average 850-^m flux 
of 0.52 ± 0.19 mJy (corresponding to a 870-/im 
flux of 0.48 ± 0.18 mJy) from a sample of 112 
-?Ceg a < 20 sBzK galaxies selected within the part 
of the Subaru/XMM-Newton deep field (SXDF) 
covered by the SCUBA HAlf Degree Extragalac- 
tic Survey (SHADES). Finally, D08 stacked the 
850- (im signal from 1421 Jf voga < 21.7 sBzK galax- 
ies and found 0.53 ± 0.06 mJy (corresponding to 
0.49 ± 0.06 mJy at 870-^m). D08 also measured 
an average 850-/j,m flux of 0.22 ± 0.18 mJy (or 
0.20 ± 0.17 mJy at 870-^ni) for a sample of 147 
-Kvega ^ 21.7 pBzK galaxies. 

For the EROs and DRGs we find stacked 870- 
/im fluxes of 0.30 ± 0.03 mJy and 0.32 ± 0.04 mJy 
respectively. Webb et al. (2004) used the Canada- 
UK Deep Submillimcter Survey 03hr and 14hr 
fields (CUDSS03 and CUDSS14, respectively) 
to perform a 850-//m stacking analysis of 164 
■FCcga S 20.7 EROs in the two fields, and found 
a stacked signal for the entire ERO sample of 
(S STOfim ) = 0.52 ± 0.09 mJy. Similarly, Takagi et 
al. (2007) measured a stacked 870-/im signal of 
0.49 ± 0.16 mJy from a sample of 201 iT voga < 20 
EROs selected within SXDF/SHADES. Turning 
to DRGs, Knudsen et al. (2005) obtained a stacked 
850- /im signal of 0.74 ± 0.24 mJy from a sample 
of 25 if voga < 22.5, submm-faint DRGs (uncor- 
rected for an average gravitational lens amplifi- 
cation of 20%). Converting to a 870-/mi flux 
density and correcting for the lens amplification 
yields 0.54 ± 0.18 mJy. Takagi et al. (2007) us- 
ing significantly shallower SCUBA maps failed 
at detecting a significant 850-/zm signal from an 
average of 67 if voga < 20, submm-faint DRGs 
((SsTo^m) = 0.39 ±0.23 mJy). 

We conclude that, within the errors, the pre- 
vious stacking studies agree well with our results. 
We also note that our study provide the first ro- 
bust (> 5-c) detection of submm-faint DRGs. 

From the stacked submm fluxes we are able to 
estimate average IR luminosities and star forma- 



4 We scale 850-fim fluxes to 870-fim fluxes assuming an op- 
tically thin, modified black body law with /3 = 1.5, such 
that: Sg7o^m /S 8 50Mm = (850/ i m/870Atm) 2 +' 3 = 0.92. 



tion rates (Table [3]) . IR luminosities are derived 
by adopting the IR-to-submm SED of Arp220 
and scaling it to the stacked submm fluxes (at 
the median rcdshifts derived from the redshift 
distributions in § I4.2p and integrating it from 
8 — 1000 /im. For comparison we also derive IR lu- 
minosities assuming that the SEDs are described 
by modified black-body law with a dust tem- 
perature of T d = 35 K and /3 = 1.5, and inte- 
grating from 8 — 1000 ^m. Star formation rates 
(SFR) are derived following Kennicutt (1998): 
SFR[M Q yr" 1 ] = 1.73 x 1O- 1O L IR [L ]. This con- 
version assumes a Salpeter initial mass function 
(Salpeter 1955). Of course, we stress that consid- 
erable uncertainty is associated with the derived 
IR luminosities and SFRs since they depend on 
the assumed SED and IMF. 

The average IR luminosities and star formation 
rates estimated here for -ftT voga < 20 sBzK, ERO 
and DRG galaxies on the basis of their stacked 
submm fluxes lie in the ranges ~l-6x 10 11 L Q 
and ~ 20 — 110 M Q yr _1 , i.e. comparable to those 
of LIRGs and ULIRGs in the local Universe. For 
the sBzK galaxies, the average IR luminosity and 
star formation rate derived here are fully consis- 
tent with UV, 24-/im and radio studies of these 
galaxies (Daddi et al. 2007). We find that ERO 
and DRG populations have significantly lower IR 
luminosities (by ~ 40%) than the sBzK galaxies. 
This is in part due to the fact that we have made 
no attempt to weed out passive EROs/DRGs in 
the stacking analysis, and the stacked submm flux 
from dusty, starforming EROs/DRGs is likely to 
be significantly higher. 

6.2. Stacking in redshift bins 

Using the photometric redshifts obtained in 
§ 14.21 we can stack our samples into separate red- 
shift bins, thereby allowing us to determine which 
redshifts are contributing the most to the stacked 
submm signals. Redshift bins were chosen so that 
they were larger than the typical redshift uncer- 
tainty, and provided roughly the same number of 
sources in each bin. The latter ensured that the 
same sensitivity was reached in each bin, thus al- 
lowing for a direct comparison. Fig. [8] shows the 
average flux densities of the different samples as a 
function of redshift. We stress that due to the es- 
sentially flat selection function of submm surveys 
over the redshift range 1 <, z <, 8 (Blain & Longair 
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Table 3 

The average IR luminosities and star formation rates of BzK, ERO and DRG galaxies, 

DERIVED FROM THEIR STACKED 870-/iM FLUXES. 



Galaxy type 


<z> L 


f R (Arp220) 


SFR b (Arp220) 


L TR < T d = 35 K ) 


SFR b (T d = 35 K) 


L IB. < T d = 30 K ) 


SFR b (T d = 30 K) 






[XIO^L© ] 


[Mq yr" 1 


[xIO^Lq ] 


[Mq yr" 1 ] 


[x 10 11 L Q ] 


[Mq /yr" 1 ] 


sBzK 


1.8 


6.3 ± 0.4 


109 ± 7 


3.7 ± 0.3 


65 ± 5 


2.1 ± 0.2 


37 ± 3 


pBzK 


1.6 


2.9 ±1.0 


50 ± 17 


1.7 ± 0.7 


30 ± 11 


1.0 ± 0.4 


17 ± 6 


ERO 


1.3 


3.6 ± 0.3 


62 ± 5 


2.3 ± 0.2 


39 ± 3 


1.3 ±0.1 


22 ± 2 


DRG 


1.4 


3-6 ± 0.4 


62 ± 7 


2.3 ± 0.3 


40 ± 4 


1.3 ±0.1 


22 ± 3 


sBzK (z < 1.4) 


1.1 


3.9 ± 0.8 


67 ± 14 


2.4 ± 0.5 


42 ± 9 


1.3 ± 0.3 


10 ± 2 


sBzK (z > 1.4) 


2.0 


5.8 ± 0.5 


101 ± 9 


3.9 ± 0.3 


67 ± 6 


2.3 ± 0.2 


39 ± 4 


pBzK (z < 1.4) 


1.2 














pBzK (z > 1.4) 


1.9 


4.7 ± 1.2 


82 ± 21 


3.0 ± 0.8 


52 ± 14 


1.7 ± 0.5 


30 ± 8 


ERO (z < 1.4) 


1.1 


2.5 ± 0.4 


43 ± 7 


1.5 ±0.2 


26 ± 4 


0.8 ±0.1 


14 ± 2 


ERO (z > 1.4) 


2.1 


5.8 ± 0.5 


101 ± 9 


3.7 ± 0.3 


64 ± 6 


2.2 ± 0.2 


37 ± 3 


DRG (z < 1.4) 


1.1 


2.8 ± 0.6 


48 ± 10 


1.7 ± 0.3 


29 ± 6 


0.9 ±0.2 


16 ± 3 


DRG (z > 1-4) 


2.4 


5.6 ± 0.6 


96 ± 10 


3.8 ± 0.4 


67 ± 7 


2.3 ± 0.2 


40 ± 4 


a IR luminositic 


s arc obtaii 


led by intcgratii 


lg trie SED over the ' 


wavelength range 8 - 


1000 rtm. 






''Star formation 




-ived using SFR[Mg yr _1 ] = 1.73 X 


10- 10 L IR [L o ] (Ke, 


micutt 1998). 







1993), the comparison of stacked submm fluxes at 
different redshifts directly translates into a com- 
parison between far-IR luminosities (and thus star 
formation rates - Kennicutt 1998) between these 
redshift bins (to the extent that submm flux is a 
measure of far-IR luminosity). 

The average submm signal from K voga < 20 
galaxies is found to be roughly constant (~ 0.1 — 
0.2 mJy) out to z ~ 1.4, and consistent with 
the average flux density of the full sample (Ta- 
ble^. At z ~ 1.7, however, the submm signal 
has increased to ~ 0.4 mJy. The stacked submm 
fluxes of z < 1.4 and z > 1.4 K vcsa < 20 sources 
are 0.15 ± 0.01 mJy and 0.47 ± 0.03 mJy, respec- 
tively (Table [2]). A similar increase in the aver- 
age submm flux at z > 1.4 is seen for the ERO 
and DRG populations. The average IR luminosi- 
ties and star formation rate for z > 1.4 EROs 
(L IR ~ 4.9 x 10 n L o and SFR ~ 85 M yr" 1 ) 
are about 3x higher than for z < 1.4 EROs 
(Lib ~ 1.8 x 10 n L Q and SFR ~ 31 M© yr" 1 ). 
For DRGs the difference is about a factor of two: 
z > 1.4 DRGs have on average L 1R ~ 4.9xlO n L 
and SFR ~ 85 M Q yr" 1 while z < 1.4 DRGs have 
L IR ~ 2.4 x 10 11 Lq and SFR ~ 42 M yr" 1 ). 
These findings also fit with the ERO and DRG 
redshift distributions (Fig. 2]), where we found ev- 
idence for two sub-populations separated at z ~ 
1.6. Due to our magnitude cut-off at if V ega < 20, 
our samples are biased towards increasingly more 
massive galaxies at higher redshifts. The above 
results, therefore, suggest that the star forma- 
tion activity in massive galaxies increase with red- 



shift, and that a dominant fraction of massive 
EROs/DRGs at z > 1.4 are actively starform- 
ing galaxies, and not old, red galaxies. The less 
massive EROs and DRGs, which aren't picked up 
at high-z given the K vesiL < 20 limit, have lower 
star formation rates, and possibly a higher frac- 
tion of evolved, passive galaxies, than the more 
massive EROs/DRGs. The gradual drop in the 
average submm signals at z > 2.5 is likely a reflec- 
tion of the incompleteness in our if-band selection 
at these redshifts. 

Turning to the BzK galaxies, we find that 
sBzKs exhibit a positive and, within the error 
bars, constant 870-/im signal (~ 0.5 mJy) over the 
redshift range 1 < z < 3. This supports the no- 
tion that the sBzK-criterion selects star forming 
galaxies across this redshift range, and that for 
our X-band magnitude limit of K vasa < 20, the 
distribution of star formation rates of sBzK galax- 
ies is roughly constant within this redshift range. 
In contrast, the pBzK galaxies show no evidence 
of significant submm signal in any redshift bin, 
which is consistent with these galaxies being de- 
void of significant star formation. D08 reported 
a 850-/mi signal of 0.89 ± 0.34 mJy (~ 2.6-cr) for 
pBzKs at z < 1.4 but no significant signal for 
pBzKs at z > 1.4. They argued that the submm 
signal for z < 1.4 pBzK galaxies was due to con- 
tamination by star forming galaxies at z < 1.4. 
In comparison, we find no significant stacked 870- 
fim signal (-0.14 ± 0.17 mJy) for z < 1.4 pBzKs 
and only a marginal signal for z > 1.4 pBzKs 
(0.36±0.12mJy). 
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Fig. 8. — Stacked 870- /urn flux densities (left panels) and the corresponding contribution to the 870-/im 
extragalactic background light (right panels) as a function of redshift for the if veg a < 20, BzK, ERO and 
DRG samples. The gray histograms in the left hand panels indicate the number of sources in each bin 
(corresponding y-axis is on the right hand side). Broader, more light-grey histograms have been used for the 
pBzK galaxies. 



6.3. Stacking in 24-/im bins 

Using the 24-^m source catalogue from the FI- 
DEL survey (which includes all sources > 27 /iJy, 
5-ct point-source sensitivity - Dickinson et al. in 
prep), a total of 466/53, 511 and 400 sources in 
the sBzK/pBzK, ERO and DRG samples, respec- 
tively, were identified at 24- /im. Their average 
submm fluxes are given in TableOalong with those 
of the 24-/xm faint (< 27/iJy) subsets. The 24-/im 
detected subsets have £ 5x higher average submm 
flux densities than the 24-/xm faint sources. This is 
not surprising since the mid-IR is known to trace 
the thermal dust emission, and one might even ex- 
pect a correlation between the submm and mid-IR 
emission. 

In order to investigate the latter, we mea- 
sured the stacked 870-/im signal as a function 
of 24- /im flux density bins (Fig. |5J). The sBzK, 
ERO and DRG galaxies all exhibit a similar be- 
haviour, namely a significant linear correlation be- 
tween the stacked 870-/im and 24- nm. flux den- 
sity up to S Mim ~ 350 fjjy, given by (S STOfim ) = 



4.5 x 10 _3 (5 24Mm ), followed by a flattening of the 
relation for 5 , 24flm > 350 /iJy. Not surprisingly, the 
pBzKs are not significantly detected at 870-/im in 
any of the 24- /im bins, and no submm-mid-IR cor- 
relation is seen. 

Although the 24-/xrn selection function (Daddi 
et al. 2007) depends strongly on redshift (unlike 
the selection function at submm wavelengths) , the 
observed S 2illln — S %70flla correlation strongly sug- 
gests that 24-/zm measurements (with S l 2 4 (jmJy Z 
350 /iJy) trace systems dominated by star forma- 
tion, and can be used to derive IR luminosities 
and star formation rates. The constant S , S70tim — 
1.3 mJy ratio for 5 , 24Mm ~ 350 — 1000 /iJy is likely 
to reflect significant contamination by an AGN at 
these high 24-/im flux densities. Since the mid- 
IR will be more sensitive to warm dust (T d ~ 
80 — 200 K) than the submm, which largely traces 
cold dust (T d ~ 20 - 60 K), the ability of the ob- 
served 24- /im emission to reliably trace IR lumi- 
nosity and star formation may be compromised 
by AGN- heated warm dust. The turn-over from 
a linear to a flat 5 870um — 5 , 24om relation occurs 
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Fig. 9. — Stacked 870-^m flux densities (left panels) and the corresponding contribution to the 870-^im 
extragalactic background light (right panels) as a function of 24-/im flux density for the BzK, ERO and 
DRG samples. The gray histograms in the left hand panels indicate the number of sources in each bin 
(corresponding y-axis is on the right hand side). Broader, more light-grey histograms have been used for the 
pBzK galaxies. An all cases, a significant S , 870flm — S , 24tim correlation is found for S 2ilim ^ 350 /iJy, indicating 
that across this flux range, 24-/im observations trace star formation. For <; 350 ^Jy, however, the 

relation flattens, suggestive of an increase in the fraction of AGN-dominated systems. 



at 5 l 24^ m — 350 /iJy, which corresponds to L m < 
1.5 x 10 12 L Q at z ~ 2. 

Our findings are in line with those of Papovich 
et al. (2007), who found that IR luminosities de- 
rived from 24-, 70-, and 160-/xm Spitzer data (of a 
sample of if-selected galaxies at 1.5 S> z 2.5 with 
S 2 4 pmjj ~ 50 — 250 /iJy) were in good agreement 
with those derived from 24-/xm data alone. For 
sources with 5 l 24)jm j y > 250 /zJy, however, the lat- 
ter would be ~ 2 — 10 x higher, suggesting that the 
AGN may contribute significantly to the high 24- 
fxm emission. Similarly, Daddi et al. (2004, 2005) 
found that a large fraction (> 30%) of sBzK galax- 
ies with L m > 1.5 x 10 12 L Q (which corresponds to 
the IR luminosity where we find a turn-over in the 
£>870vm — S , 24 M m relation) show an excess of emission 
in the near-IR (rest-frame) and are statistically de- 
tected in hard X-rays - evidence of powerful AGN 
in these very IR-luminous systems. 



6.4. Contributions to the extragalactic 
background light 

Turning to the contribution to the submm ex- 
tragalactic background light (EBL) by the differ- 
ent populations, we adopt the spectral approxi- 
mations to the EBL at submm wavelengths from 
COBE/FIRAS (Puget et al. 1996; Fixsen et al. 
1998), including their uncertainties. In doing 
so we adopt a value of the EBL at 870-/im of 
44 ± 15 Jydeg~ 2 . From the surface densities of the 
samples we derive their contributions to the 870- 
\im. extragalactic background light (Tabled]). We 
find that the total contribution from all -fsT vcga < 20 
sources to the 870-^m EBL is 6.61 ± 0.34 Jydeg" 2 , 
or 15.0 ± 5.2 %. For the BzK galaxies we find 
that the passive ones contribute < 1 % to the 
870- /im EBL, while the star forming population 
contribute 1.70 ± 0.24 Jydeg -2 , corresponding to 
3.9±1.4 %. The EROs contribute a similar amount 
to the EBL at 870-/xm as the sBzK galaxies (1.85± 
0.20 Jydeg -2 or 4.2 ± 1.5%), while the contribu- 
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tion from DRGs, owing to their lower surface den- 
sity, is about 80% smaller (1.09 ± 0.18 Jydeg" 2 or 
2.5±0.9 %). The combined BzK/ERO/DRG sam- 
ple, which includes the bulk of massive starforming 
galaxies at z 'z 1, contribute with 8.4 ± 2.9% to 
the EBL at 870-/im after accounting for overlap 
between the populations. 

According to Takagi et al. (2007) sBzK galax- 
ies down to A vcga ^ 20 contribute with 3.8 ± 
1.2 Jydeg" 2 (or 8.3 ± 3.9%) to the background 
light at 850-^m (46 ± 16 Jydeg" 2 - Puget et al. 
1996; Fixsen et al. 1998). Webb et al. (2004) found 
that the ERO population down to A" vcga < 20.7 
constitutes 7-11% of the EBL at 850- /im, while 
Takagi et al. (2007) reported a 850- /im EBL contri- 
bution from A voga £ 20 EROs of 5.1 ± 1.5 Jy deg~ 2 
(11 ±5%) to the total background at 850-/Ltm. 
Knudsen et al. (2005) found that DRG galaxies 
down to A voga < 22.5 contribute ~ 7.7Jydeg -2 
(~ 17 %) of the EBL at 850-^m. 

Next, let us look at the contribution to the EBL 
at 870-/im from the different samples as a function 
of redshift and 24- /im flux bins (right panels in 
Fig. [5] & [9]). We calculate the contributions using 
the stacked submm flux and the surface density of 
sources in each redshift and 24-/im bin. 

Looking at the redshift dependence first, we see 
that the strongest contribution to the 870-^m EBL 
by A vcga < 20 galaxies in a given redshift bin, 
is coming from sources in the redshift range 1 < 
z < 2 - they contribute with 7.3 ± 0.2 Jydeg -2 
to the submm EBL, which corresponds to ~ 50% 
of the total contribution from AT voga < 20 galaxies 
(Table O. Although A voga < 20 galaxies at z > 
2 have significant submm emission (§ 16. 2p . their 
low abundance means that they contribute ^ 1% 
to the observed submm EBL. In contrast, z < 1 
AT vega < 20 galaxies are abundant (Fig. [4j, but 
their low average submm fluxes implies a < 1 % 
contribution to the submm EBL. 

The bulk (~ 80 %) of the contribution to the 
submm EBL from 24-fim detected (i.e. S 24l ,., y > 
27 /iJy) sBzK, ERO and DRG galaxies comes from 
sources with S 24llJy ~ 50 — 350 /iJy, with only a 
minor fraction coming from brighter, presumably, 
AGN-dominated sources. 



6.5. Stacking across the BzK and RJK di- 
agrams 

The sBzK and pBzK selections are designed to 
locate star forming and passive galaxies in the 
redshift range 1.4 < z < 2.5, while the ERO 
and DRG colour criteria are designed to select ex- 
tremely red (either due to dust extinction or old 
age) galaxies at z > 1. Still, the exact defini- 
tions of these colour criteria are somewhat arbi- 
trary. Given the large number of sources avail- 
able to us, we are in a position to construct stacks 
of the 870-/j,m signal from statistically significant 
subsets of galaxies across the BzK and RJK di- 
agrams, thus allowing us to see where in these 
colour-colour diagrams the submm signal is com- 
ing from. Ultimately, this may allow us to fine- 
tune the sBzK/pBzK criteria, as well as poten- 
tially identify regions of the RJK-diagram contain- 
ing dusty/starforming vs. old/passive EROs and 
DRGs. 

In Fig . [TOl and [Til we show the signal- and S/N- 
contours of the stacked 870-^m signal obtained 
across the BzK- and RJK-diagrams, where the 
stacking has been carried out within equally sized 
grid-cells. We have done this for the full sam- 
ple, as well as for subsets of the sample within 
the redshift intervals < z < 0.8, 0.8 < z < 1.5, 
and 1.5 < z < 3.0. These intervals were moti- 
vated by the redshift distributions in Fig.|4j which 
showed evidence of distinct populations between 
0.8 < z < 1.5 and 1.5 < z < 3.0. 

Considering Fig. [T0k first, we see that the sBzK 
criterion by Daddi et al. (2004) seems to be robust 
in the sense that the bulk of the 870-/im signal is 
emerging from the sBzK region. Comparing with 
Fig. [1^, indicates that virtually the entire submm 
signal is coming from the subset of sBzK galaxies 
(with a slight overlap into the pBzK and non-BzK 
regions), which have also been classified as EROs. 
Due to their extremely red colours and relatively 
strong submm signal, these galaxies are very likely 
to be amongst the most dusty, star forming sources 
of the submm-faint A vcga < 20 selected galaxies, 
yet their brightness ensures that they are detected 
in the blue (thus qualifying as sBzK galaxies). 

The submm signal clearly extends into the 
pBzK region suggesting that some contamination 
by star forming galaxies occurs. This is partly 
due the low density of sources in the pBzK-region 
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Fig. 10. — Contours of the stacked 870-/im signal (in units of mJy) and S/N across the BzK diagram, 
obtained by stacking the submm flux in a regular grid (shown in blue). We have stacked the entire sample 
within the full redshift range, as well as subsets within the redshift ranges z = — 0.8, 0.8 — 1.5, and 
1.5 — 3.0. The sBzK- and pBzK-criteria by Daddi et al. (2004) are shown as solid green lines, while the 
refined selection criteria proposed in this paper are shown in red. The contours clearly supports the BzK- 
selection technique as the submm signal is seen to be almost entirely dominated by sBzK galaxies in the 
redshift range z — 1.5 — 3. 



which implies that more sources will scatter into 
the region (due to photometric errors) than out of 
it. Based on the 870-/im signal and S/N contours 
in Fig. [ini we propose a refinement of the sBzK 
and pBzK criteria, namely 

{z-K) AB > (B-z) AB - 0.5 (6) 
(B-z) AB < 3.4 (7) 



for sBzKs, and 



(z-K) AB > 2.9 
(B-z) AB > 3.4 



(8) 
(9) 



for pBzKs. 

That these refined selection criteria might do 
a better job at selecting starforming and passive 



BzK galaxies, is confirmed if we consider how the 
submm signal across the BzK-diagram changes 
with redshift (Fig. [TOb-h). We find that only 
marginal 870-/im emission emerges from galaxies 
in the redshift interval < z < 0.8, while in the 
0.8 < z < 1.5 interval, significant submm emis- 
sion (~ 0.6 mJy) starts to appear from galaxies 
lying within the sBzK region. At 1.5 < z < 3 
a strong submm signal (~ 1 mJy) and nearly all 
of the submm emission is coming from the sBzK 
region. This observed migration of the stacked 
submm signal towards the sBzK region as a func- 
tion of redshift, illustrates the ability of the sBzK- 
criterion to select starforming galaxies at z > 1.5. 

Turning to the RJK-diagram, the bulk of the 
submm-signal is found to come from the overlap 
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Fig. 11. — Contours of the stacked 870-^m signal (in units of mJy) and S/N across the RJK diagram, 
obtained by stacking the submm flux in a regular grid (shown in blue). We have stacked the entire sample 
within the full rcdshift range, as well as subsets within the redshift ranges z = — 0.8, 0.8 — 1.5, and 1.5 — 3.0. 
The ERO- and DRG-criteria are shown as vertical and horizontal solid green lines, while the starforming vs. 
passive ERO criterion proposed by Pozzetti & Mannucci (2000) is illustrated by the green dashed line (with 
starforming and passive EROs lying above and below the line, respectively). The strongest submm signal is 
coming from galaxies in the redshift range z = 1.5 — 3 lying in the DRG-ERO overlap region. 



region between DRGs and EROs, in particular 
from ERO/DRG galaxies with (J - K) AB > 1.9. 
The latter, together with the tentative evidence of 
an underdensity of sources at {J — K) AB — 1.9 (see 
Fig. [1]), suggests that (J — K) AB > 1.9 may be a 
useful criterion for selecting a dusty, starforming 
galaxies. This is strengthened by the fact that the 
submm signal becomes even stronger and further 
concentrated towards the ERO/DRG overlap re- 
gion if we consider only the subset of sources at 
z = 1.5 — 3 (Fig. [TTk-f0. In the redshift range 
0.8 < z < 1.5, the submm signal is much weaker 
(~ 0.6 mJy) and comes from ERO/DRG galaxies 
with (J-K) AB < 1.9. 

It is well-known that the i\~ voga < 20 ERO pop- 
ulation is a heterogeneous population, consisting 



of roughly a fifty-fifty mix of dusty, star forming 
EROs and evolved, passive EROs (Spinrad et al. 
1999; Dey et al. 1999; Cimatti et al. 1999; Man- 
nucci et al. 2002). Pozzetti & Mannucci (2000) 
argued that a crude separation between the two 
types of EROs could be made based on the cri- 
terion (J - K) AB = 0.M(R - K) AB - 0.22 (and 
(R — Ks)ab > 3.35 jH, with dusty, star forming 
EROs lying above the relation and old, passive 
EROs below it. While our stacking analysis lends 
some merit to the Pozzetti & Mannucci criterion, 
as the bulk of the submm signal is clearly found 
above it, significant submm emission is also de- 
tected in the passive ERO region (in particular for 



5 We have converted the Pozzetti &; Mannucci criterion, 
which was in the vega system, into AB magnitudes. 
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sources at z > 1.5), suggesting that blindly apply- 
ing the criterion does not produce clean samples 
of starforming and passive EROs (see also Smail 
ct al. 2002b). 

Finally we caution that by adopting equally 
sized grid cells, some cells will contain a sig- 
nificantly larger number of sources than others, 
thereby introducing a potential skewing of the 
measured S/N across the diagrams (which is why 
we also show the variation of the average submm 
flux density across the diagrams). As a check on 
our results, we therefore adopted two alternative 
methods for binning the sources. First, an adap- 
tive mesh was applied by requiring that no more 
than 200 sources were allowed within a given cell, 
and, secondly, the 50 nearest neighbours of each 
source were identified and stacked. Reassuringly, 
the two additional binning methods gave results 
in good agreement with the regular grid results. 

7. Summary 

Using the APEX/LABOCA 870-^m map of 
the ECDF-S (Weifi et al. 2009) along with the 
publically available MUSYC survey data of this 
field (Taylor et al. 2008), we have performed a 
submm stacking analysis of 8266 if-band selected 
(-ft'vega < 20) galaxies, as well as subsets of 737 
DRG, 1253 ERO, and 744/149 sBzK/pBzK galax- 
ies. Photometric rcdshifts have been derived for 
the full if voga < 20 sample using UBVRIzJHK 
data from MUSYC, thereby allowing us to study 
stacked submm flux densities as a function of red- 
shift. This represents the largest submm stacking 
analyses of near-IR selected galaxies to date. The 
main results are summarized below: 

• We measure stacked 870-^m signals of 
0.20 ± 0.01 mJy (20.0cr), 0.45±0.04mJy (11.3cr), 
0.42 ± 0.03 mJy (U.Oa), and 0.41 ± 0.04mJy 
(10.3<t) for the X vcga < 20, BzK, ERO and 
DRG samples, respectively. Splitting the BzK 
galaxies up into star-forming (sBzK) and passive 
(pBzK) galaxies, the former is significantly de- 
tected (0.48 ± 0.04 mJy, 12.0a) while the latter, as 
expected, is not (0.27 ± 0.10 mJy, 2.7a). This im- 
plies that i"T vcga < 20 galaxies are responsible for 
15 ± 5% of the EBL at 870-^m. sBzK galaxies, 
EROs and DRGs (brighter than K vosa < 20) are 
found to contribute ~ 4% of the 870-^m back- 



ground each. 

• Performing the stacking analysis in redshift 
bins, it is found that the stacked submm signal 
from if vcga < 20 galaxies, as well as the ERO and 
DRG sub-samples, is coming from sources in the 
redshift range 1.4 < z < 2.5, while for BzK galax- 
ies the signal remains constant with redshift. As- 
sessing the contribution to the submm EBL from 
the different samples function of redshift is 
complicated by the fact that these are flux-limited 
samples, although, we can conclude that the bulk 
of the submm light from i^ vega < 20 galaxies comes 
from sources at z > 1. 

• We find a linear correlation between stacked 
submm flux and 24-/xm flux density for sBzK, 
ERO and DRG galaxies with ^ 350 /zJy. 
This correlation suggests that the 24-^m emis- 
sion from Si^m <, 350 /iJy galaxies is domi- 
nated by star formation, and consequently can 
be used as a robust tracer of star formation. At 
<S l 24/im > 350 /xJy we find that the stacked 870-/xm 
flux density becomes constant (~ l.SmJy) and 
independent of S^m, which is likely due to AGN 
starting to contribute significantly to the 24-/im 
as well as 870-/im emission. 

• In an effort to isolate a subset of BzK, ERO, 
and DRG galaxies responsible for the bulk of the 
stacked 870-^m emission we have measured the 
significance of the stacked submm signal across 
the BzK and RJK diagrams, identifying the re- 
gions of strongest submm emission. We find that 
the stacked submm signal from submm-faint sBzK 
galaxies is dominated by the subset of sources 
which also fulfill the ERO criterion. These are 
likely to be dusty, starforming galaxies, which are 
sufficiently bright in the blue to be selected as 
sBzK galaxies. The majority of these sources are 
found in the redshift range z = 1.5 — 3, in line with 
the BzK-selection criterion proposed by Daddi et 
al. (2004) . Guided by the stacked submm-contours 
we propose slightly modified BzK-selection crite- 
ria, namely 

(z-K) AB > (B-z) AB -0.5 (10) 
(B-z) AB < 3.4 (11) 

(12) 
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for sBzKs, and 

(z-K) AB > 2.9 (13) 
(B-z) AB > 3.4 (14) 

(15) 

for pBzKs. 

In the RJK-diagram we find that the strongest 
submm signal comes from galaxies in the ERO /DRG 
overlap region with (J — K) AB > 1.9, which are 
predominantly found at z > 1.5. 
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